/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/* Direct evaluation of the radial distribution function. Adapted from [Rapaport, p.90].
   Sets of diffusion measurements are begun at fixed time intervals, so that several 
   sets of measurements based on different time origins will be in progress simultaneously
   because of the overlapped measurements.
   Memory is allocated in module startUp.c, function AllocArrays(). */

#include "dataStructs.h"

extern const RMat *siteRotMat;
extern const VecR region, *siteCoords; 
extern const Site *site; 
extern const int nSites, limitWatRdf, sizeHistWatRdf, stepWatRdf, rdfBallsOnly;
extern const real deltaT, rangeWatRdf, timeNow;
extern real *histWatRdfOO, *histWatRdfOH, *histWatRdfHH;
extern int countWatRdf;

void ZeroWatRdfHistograms(); /* external function declaration */

static void PrintWatRdf()
{
  real distance; // atom-atom distance -- the x-axis of final plot
  int n; // counter

  FILE *rdfFilePtr;
  char rdfSnapFile[ 35 ];
  char time_ps[ 9 ]; // string containing timeNow, e.g. "100", in ps

  sprintf( time_ps, "%d", ( int ) ( 0.5 + timeNow * TIME_ps ) );
  rdfSnapFile[ 0 ] = '\0';
  strcat( rdfSnapFile, "waterRdf_" );
  strcat( rdfSnapFile, time_ps );
  strcat( rdfSnapFile, "ps" );
  strcat( rdfSnapFile, ".dat" );
  rdfFilePtr = fopen( rdfSnapFile, "w" );
  
  for ( n = 0; n < sizeHistWatRdf; n++ ) {
    distance = ( n + 0.5 ) * rangeWatRdf / sizeHistWatRdf;
    fprintf( rdfFilePtr, "%8.3f", distance );
    fprintf( rdfFilePtr, "%10.6f", histWatRdfOO[ n ] );
    fprintf( rdfFilePtr, "%10.6f", histWatRdfOH[ n ]  );
    fprintf( rdfFilePtr, "%10.6f", histWatRdfHH[ n ]  );
    fprintf( rdfFilePtr, "\n" );
  }
  fflush( rdfFilePtr );
}

void EvalWatRdf() // [Rapaport, p 223]
{
  VecR dr, shift, r;
  real deltaR, normFac, rr, rangeRdfSquared;
  int i, j, n;
  RMat rMat;
  VecR OH1b, OH2b, /* oxygen-hydrogen vectors in body-fixed frame */
    OOb, t, 
    H1is, H2is, Ois, // resolved position of hydrogens + oxygen in the spaced-fixed frame
    H1js, H2js, Ojs; // resolved position of hydrogens + oxygen in the spaced-fixed frame
  
  // water internal geometry in Brahms units
  VSet( OH1b, 0.,  0.075 / LENGTH_nm, 0.053 / LENGTH_nm ); /* from Bratko et al, J Chem Phys, 83, 6367, 1985 */
  VSet( OH2b, 0., -0.075 / LENGTH_nm, 0.053 / LENGTH_nm );
  VSet( OOb,  0.,  0.,             -0.00654 / LENGTH_nm );

  deltaR = rangeWatRdf / sizeHistWatRdf;

  rangeRdfSquared =  Sqr( rangeWatRdf );

  for ( i = 0; i < nSites - 1; i++ ) {
    for ( j = i + 1; j < nSites; j++ ) {
      VSub( dr, siteCoords[ i ], siteCoords[ j ] );
      VZero( shift ); // zero auxiliary vector "shift"
      VShiftAll( dr ); // if MIC.x, shift = \pm region.x, etc.
      VWrapAll( dr );
      rr = VLenSq( dr );
      if ( rr < rangeRdfSquared ) {

	MCopy( rMat.u, siteRotMat[ i ].u );      
	VCopy( r, siteCoords[ i ] );
	MVMulT( t, rMat.u, OH1b );
	VAdd( H1is, r, t );
	MVMulT( t, rMat.u, OH2b );
	VAdd( H2is, r, t );
	MVMulT( t, rMat.u, OOb );
	VAdd( Ois, r, t );

	MCopy( rMat.u, siteRotMat[ j ].u );      
	VCopy( r, siteCoords[ j ] );
	MVMulT( t, rMat.u, OH1b );
	VAdd( H1js, r, t );
	MVMulT( t, rMat.u, OH2b );
	VAdd( H2js, r, t );
	MVMulT( t, rMat.u, OOb );
	VAdd( Ojs, r, t );

	// computing Oxygen-Oxygen rdf
	VSub( dr, Ois, Ojs);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfOO[ n ];
	}
	
	// computing Oxygen-Hydrogen rdf
	VSub( dr, Ois, H1js);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfOH[ n ];
	}
	VSub( dr, Ois, H2js);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfOH[ n ];
	}
	VSub( dr, Ojs, H1is);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfOH[ n ];
	}
	VSub( dr, Ojs, H2is);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfOH[ n ];
	}

	// computing i-j Hydrogen-Hydrogen rdf
	VSub( dr, H1is, H1js);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfHH[ n ];
	}
	VSub( dr, H1is, H2js);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfHH[ n ];
	}
	VSub( dr, H1js, H2is);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfHH[ n ];
	}
	VSub( dr, H2js, H2is);
	VWrapAll( dr );
	rr = VLenSq( dr );
	if ( rr < rangeRdfSquared ) {
	  n = sqrt( rr ) / deltaR;
	  ++histWatRdfHH[ n ];
	}
      }
    }
  }
  
  ++countWatRdf;
  if ( countWatRdf == limitWatRdf ) {
    normFac = VProd( region ) / ( 2. * PI * Cube( deltaR ) * Sqr( nSites ) * countWatRdf );
    for ( n = 0; n < sizeHistWatRdf; n++ ) {
      histWatRdfOO[ n ] *= normFac / Sqr( n - 0.5 );
      histWatRdfOH[ n ] *= .25 * normFac / Sqr( n - 0.5 );   // (/ 4.) because of the 4 H-H pair for every O-O pair [Rapaport, p.224]
      histWatRdfHH[ n ] *= .25 * normFac / Sqr( n - 0.5 );   //  ditto
    }
    PrintWatRdf();
    countWatRdf = 0;
    ZeroWatRdfHistograms();
  } 
}
